Patients with severe asthma can be treated with biologic medicinces. NICE has approved four biologics for use in England for people with severe asthma; benralizumab, mepolizumab and most recently reslizumab for treating severe eosinophilic asthma and omalizumab for treating severe persistent confirmed allergic IgE‑mediated asthma. Our team delivers OpenPrescribing.net, a publicly funded and openly accessible explorer for NHS primary care prescribing supported by prescribing data openly published by the NHS Business Services Authority. We were concerned that the NHS did not share hospital medicines data in a similar manner to primary care and that it should be shared (see Goldacre and MacKenna 2020). Asthma UK have similarly advocated that hospital medicines data should be shared to enable transparent assessment of access to new treatments and accountability for public spending.
In September 2020 NHS BSA published hospital medicines data for the first time. We have prepared the following notebook for investigating the use of asthma biologic medicines in English hospitals. If you have any feedback or insight The DataLab team can be contacted at ebmdatalab@phc.ox.ac.uk.
The list snomed_asthma_mabs_list contains the SNOMED codes for the medicines we are investigating.
# Define vector with selected SNOMED codes for Asthma Medicines
# Consider getting these from opencodelists at some point maybe
snomed_asthma_mabs_list <- c("18671311000001102",
"18671411000001109",
"37564311000001105",
"37564411000001103",
"31210311000001108",
"33999611000001101",
"34812511000001102",
"37854611000001108",
"35298511000001102")Before we can analyse trends and variation in the use of asthma biologics we need to prepare our dataset. This analysis uses three differen sources of data.
The NHS Digital “GP mapping file” which provides STP and region names mapped to STP and region ODS codes, published here.
The NHS Digital “Etr” file that maps trust organisation codes to trust names, STP ODS codes and region ODS codes, published here.
# Load ETR data
df_etr <- readr::read_csv(here::here("data/etr_tidy.csv"))
# Load stp to regions data
stp_to_region_map <- read_csv(here::here("data/gp-reg-pat-prac-map.csv")) %>%
group_by(STP_CODE, STP_NAME) %>%
summarise(COMM_REGION_NAME = first(COMM_REGION_NAME),
COMM_REGION_CODE = first(COMM_REGION_CODE)) %>%
janitor::clean_names()# Sustainability and Transformation Partnerships (STPs)
df_etr %>%
left_join(stp_to_region_map) %>%
select(ods_name, stp_name, comm_region_name) %>%
mutate(stp_name = fct_explicit_na(stp_name),
comm_region_name = fct_explicit_na(comm_region_name)) %>%
reactable::reactable(filterable = TRUE,
columns = list(ods_name = reactable::colDef(name = "Name",
minWidth = 200),
stp_name = reactable::colDef(name = "STP",
minWidth = 150),
comm_region_name = reactable::colDef(name = "Region",
minWidth = 70)
),
style = list(fontSize = "12px"),
highlight = TRUE)# Secondary Care Medicines Data
# Connect to database, filter, and collect data.
# Get SCMD dataset
# Check negative quantities, seems to be a small (0.7%) problem in the data
# dplyr::tbl(conn_ebm_scmd, "scmd") %>%
# filter(between(year_month, "2019-01-01", "2019-12-31")) %>%
# # arrange(desc(year_month))
# filter(total_quanity_in_vmp_unit > 0) %>%
# count()
# group_by(vmp_snomed_code) %>%
# summarise(n = n(),
# sum = sum(total_quanity_in_vmp_unit < 0))
db_scmd <- dplyr::tbl(conn_ebm_scmd, sql_query_scmd)
# Create dataframe for table
db_scmd_asthma <- db_scmd %>%
dplyr::filter(vmp_snomed_code %in% snomed_asthma_mabs_list)
df_scmd_asthma <- dplyr::collect(db_scmd_asthma)
# df_scmd_asthma %>% skimr::skim()
# Tidy tidy tidy data
df_scmd_asthma_names <- df_scmd_asthma %>%
dplyr::left_join(dplyr::select(df_etr, ods_code, ods_name, stp_code), by = "ods_code") %>%
# some data cleaning as scmd uses some ods codes that are not up to date
mutate(stp_code = as.character(stp_code),
stp_code = case_when(
ods_code == "RQ6" ~ "QYG",
ods_code %in% c("RNL", "RE9", "RLN") ~ "QHM",
ods_code %in% c("RM2", "RW3") ~ "QOP",
ods_code == "RGQ" ~ "QJG",
ods_code == "RJF" ~ "QJ2",
ods_code == "RR1" ~ "QHL",
TRUE ~ stp_code
))# Fill explicit missing and create dataset for sparkline in table
df_tab_sparkline <- df_scmd_asthma_names %>%
select(-c(vmp_product_name, ods_name, stp_code, stp_code, ods_name)) %>%
arrange(ods_code, vmp_snomed_code, year_month) %>%
as_tsibble(key = c(ods_code, vmp_snomed_code), index = year_month) %>%
fill_gaps(total_quantity = 0, .full = TRUE) %>%
tidyr::fill(.direction = "down") %>%
as_tibble() %>%
mutate(year_month = floor_date(year_month, unit = "month")) %>%
group_by(year_month, ods_code, vmp_snomed_code) %>%
arrange(ods_code, vmp_snomed_code, year_month) %>%
mutate(total_quantity = sum(total_quantity)) %>%
arrange(ods_code, vmp_snomed_code, year_month) %>%
distinct() %>%
group_by(ods_code, vmp_snomed_code) %>%
dplyr::summarise(count_sparkline = list(total_quantity)) %>%
group_by(ods_code, vmp_snomed_code) %>%
dplyr::mutate(total_quantity = sum(unlist(count_sparkline)))
# Create lookup datasets for joining
# SNOMED
vmp_snomed_names_lookup <- df_scmd_asthma_names %>%
select(vmp_snomed_code, vmp_product_name) %>%
dplyr::distinct()
# Trust
trust_names_lookup <- df_scmd_asthma_names %>%
select(ods_code, ods_name, stp_code) %>%
dplyr::distinct()
# Join data
df_tab_sparkline <- df_tab_sparkline %>%
ungroup() %>%
arrange(ods_code) %>%
left_join(trust_names_lookup, by = c("ods_code")) %>%
left_join(vmp_snomed_names_lookup, by = c("vmp_snomed_code")) %>%
mutate(count_box = count_sparkline)
# See the ?tippy documentation to learn how to customize tooltips
with_tooltip <- function(value, tooltip, ...) {
div(style = "text-decoration: underline; text-decoration-style: dotted; cursor: help",
tippy(value, tooltip, ...))
}
# Create table
df_tab_sparkline %>%
select(ods_name, vmp_product_name,vmp_snomed_code, count_sparkline, count_box, total_quantity,
-ods_code, -stp_code) %>%
reactable(filterable = TRUE,
defaultSorted = c("ods_name", "total_quantity"),
groupBy = c("ods_name"),
columns = list(
ods_name = reactable::colDef(name = "Trust",
# header = with_tooltip("Trust", "ods_name"),
minWidth = 200),
count_sparkline = colDef(name = "Trend",
header = with_tooltip("Trend", "Note that the y axis cannot be compared across different entries."),
minWidth = 50,
cell = function(value, index) {
sparkline(df_tab_sparkline$count_sparkline[[index]])
}),
# count_box = colDef(name = "Boxplot",
# header = with_tooltip("Boxplot", "Boxplot of monthly quantity for entries with more than 10 data points."),
# minWidth = 40,
# cell = function(value, index) {
# if (sum(df_tab_sparkline$count_box[[index]] > 0) > 10) {
#
# sparkline(df_tab_sparkline$count_box[[index]],
# type = "box")
# }
# }),
count_box = reactable::colDef(show = FALSE),
total_quantity = reactable::colDef(name = "Quantity",
# header = with_tooltip("Total quantity", "total_quantity"),
minWidth = 50,
aggregate = "sum",
format = reactable::colFormat(digits = 0)),
vmp_product_name = reactable::colDef(name = "Product",
# header = with_tooltip("Product", "total_quantity"),
minWidth = 150,
cell = function(value, index) {
vmp_snomed_code <- paste0("SNOMED: ", df_tab_sparkline$vmp_snomed_code[index])
vmp_snomed_code <- if (!is.na(vmp_snomed_code)) vmp_snomed_code else "Unknown"
div(
div(style = list(fontWeight = 600), value),
div(style = list(fontSize = 10), vmp_snomed_code))
}
),
vmp_snomed_code = reactable::colDef(show = FALSE)
),
style = list(fontSize = "12px"),
highlight = TRUE
)Information from the dm+d (Dictionary of Medicines and Devices) on the strength and vmp quantity of each asthma biologic at VMP level, using data hosted on the DataLab BigQuery server.
db_dmd_info <- dplyr::tbl(conn_ebm_scmd, sql_query_dmd_info)
df_dmd_info <- db_dmd_info %>%
filter(vmp_snomed_code %in% snomed_asthma_mabs_list) %>%
collect()We use information on the daily defined dose (DDD) of each asthma biologic, so that the volume of each VMP can be compared directly once converted to DDDs. The WHO publish DDDs online:
# Define tibble with mg_per_ddd for join later
# Is there a way this step could be automated?
lookup_mg_per_ddd <- tibble::tribble(~vtmnm, ~mg_per_ddd,
"Omalizumab", 16,
"Mepolizumab", 3.6,
"Reslizumab", 7.1,
"Benralizumab", 0.54)
df_scmd_asthma_mg <- df_scmd_asthma_names %>%
left_join(df_dmd_info, by = c("vmp_snomed_code", "vmp_product_name")) %>%
left_join(lookup_mg_per_ddd, by = "vtmnm") df_scmd_asthma_mg %>%
select(vmp_snomed_code, vtmnm, form_descr, udfs, udfs_descr,
strnt_nmrtr_val, strnt_nmrtr_uom, strnt_dnmtr_val,strnt_dnmtr_descr, mg_per_ddd) %>%
distinct() %>%
reactable(filterable = TRUE,
columns = list(
vmp_snomed_code = reactable::colDef(name = "SNOMED",
minWidth = 100),
vtmnm = reactable::colDef(name = "Name",
minWidth = 100),
form_descr = reactable::colDef(name = "Form",
minWidth = 80),
# udfs is the VMP unit dose form strength
udfs = reactable::colDef(name = "Value",
minWidth = 40),
udfs_descr = reactable::colDef(name = "Unit",
minWidth = 40),
# strnt_nmrtr
strnt_nmrtr_val = reactable::colDef(name = "Numerator",
minWidth = 60,
format = colFormat(suffix = " mg")),
strnt_nmrtr_uom = reactable::colDef(show = FALSE),
# strnt_dnmtr
strnt_dnmtr_val = reactable::colDef(name = "Denominator",
minWidth = 60,
format = colFormat(suffix = " ml")),
strnt_dnmtr_descr = reactable::colDef(show = FALSE),
mg_per_ddd = reactable::colDef(name = "mg/ddd",
minWidth = 50)),
columnGroups = list(
colGroup(name = "UDFS", columns = c("udfs", "udfs_descr")),
colGroup(name = "Strength", columns = c("strnt_nmrtr_val", "strnt_dnmtr_val"))
),
style = list(fontSize = "12px"),
highlight = TRUE
)The final data cleaning step is to convert the volume from VMP quantity (as provided in the SCMD dataset) to volume in DDDs.
df_scmd_asthma_ddd <- df_scmd_asthma_mg %>%
mutate(volume_singles = total_quantity / udfs,
volume_mg_strength = volume_singles * if_else(is.na(strnt_dnmtr_val),
true = strnt_nmrtr_val,
false = strnt_nmrtr_val *
(udfs / strnt_dnmtr_val)),
volume_ddd = volume_mg_strength / mg_per_ddd)temp_ggplot <- df_scmd_asthma_ddd %>%
group_by(year_month, vtmnm) %>%
summarise(volume_ddd = sum(volume_ddd)) %>%
ggplot(aes(x = year_month,
y = volume_ddd,
colour = vtmnm,
group = vtmnm)) +
geom_line(size = 1, alpha = 0.5) +
geom_point(aes(text = paste0("<b>Month:</b> ",
lubridate::month(year_month, label = TRUE), " ",
lubridate::year(year_month), "<br>",
"<b>Volume:</b> ", round(volume_ddd, 0), "<br>",
"<b>Medication:</b> ", vtmnm)), size = 2) +
scale_x_date(date_breaks = "4 month", date_labels = "%b %y") + scale_colour_viridis_d() +
labs(x = NULL, y = "Defined Daily Dose",
colour = NULL) +
geom_vline(xintercept = as.numeric(as.Date("2020-03-31")),
color = "orange",
linetype = 2,
lwd = .5,
alpha = .5) +
theme(text = element_text(size = 12))
# temp_ggplot
plotly::ggplotly(temp_ggplot,
tooltip = "text") %>%
plotly::config(displayModeBar = FALSE)Figure. National prescribing of asthma biologics over time.
df_scmd_asthma_ddd_map <- df_scmd_asthma_ddd %>%
left_join(stp_to_region_map, by = "stp_code")temp_ggplot <- df_scmd_asthma_ddd_map %>%
group_by(comm_region_name, year_month, vtmnm) %>%
summarise(volume_ddd = sum(volume_ddd)) %>%
mutate(comm_region_name = fct_explicit_na(comm_region_name)) %>%
ggplot(aes(x = year_month, y = volume_ddd,
colour = comm_region_name, group = comm_region_name)) +
geom_vline(xintercept = as.numeric(as.Date("2020-03-31")),
color = "orange",
linetype = 2,
lwd = .5,
alpha = .5) +
geom_line(size = 1, alpha = 0.5) +
geom_point(aes(text = paste0("<b>Month:</b> ",
lubridate::month(year_month, label = TRUE), " ",
lubridate::year(year_month), "<br>",
"<b>Region:</b> ", comm_region_name, "<br>",
"<b>Volume:</b> ", round(volume_ddd, 0), "<br>",
"<b>Medication:</b> ", vtmnm)),
size = 2) +
scale_x_date(date_breaks = "4 month", date_labels = "%b %y") +
scale_colour_viridis_d(end = 1) +
labs(x = NULL,
y = "Defined Daily Dose",
colour = NULL) +
facet_wrap(~vtmnm, ncol = 1) +
theme(text = element_text(size = 12))
# temp_ggplot
plotly::ggplotly(temp_ggplot,
tooltip = "text") %>%
plotly::config(displayModeBar = FALSE)Figure. Regional prescribing of asthma biologics over time.
temp_ggplot <- df_scmd_asthma_ddd_map %>%
select(year_month, ods_code, vtmnm, volume_ddd, comm_region_name) %>%
mutate(comm_region_name = fct_explicit_na(comm_region_name)) %>%
group_by(comm_region_name, vtmnm) %>%
summarise(volume_ddd = sum(volume_ddd)) %>%
group_by(comm_region_name) %>%
mutate(prop_use = volume_ddd / sum(volume_ddd),
pos = cumsum(volume_ddd) - volume_ddd/2,
total = sum(volume_ddd)) %>%
ungroup() %>%
mutate(comm_region_name = fct_reorder(comm_region_name, total)) %>%
ggplot(aes(comm_region_name)) +
geom_bar(aes(y = volume_ddd,
fill = vtmnm,
text = paste0("<b>Region:</b> ", comm_region_name, "<br>",
"<b>Total volume in ddd:</b> ", round(total, 0), "<br>",
# "<b>Medication:</b> ", vtmnm, "<br>",
"<b>", vtmnm , " volume in ddd (%):</b> ", round(volume_ddd, 0), " (", scales::percent(prop_use, accuracy = 0.1), ")"
)
),
stat='identity',
# position = position_dodge()
) +
scale_fill_viridis_d() +
labs(subtitle = paste0("From: ", min(df_scmd_asthma_ddd_map$year_month), " to ", max(df_scmd_asthma_ddd_map$year_month)),
x = NULL,
y = "Defined Daily Dose",
fill = NULL) +
scale_y_continuous(labels = scales::comma) +
theme(text = element_text(size = 12)) +
coord_flip()
# temp_ggplot
plotly::ggplotly(temp_ggplot,
tooltip = "text") %>%
plotly::config(displayModeBar = FALSE)Figure. Total regional prescribing of asthma biologics.
df_scmd_asthma_ddd_map_temp <- df_scmd_asthma_ddd_map %>%
filter(year_month >= as.Date("2019-08-01") & year_month <= as.Date("2020-07-01")) %>%
select(year_month, ods_code, vtmnm, volume_ddd, stp_name) %>%
mutate(stp_name = fct_explicit_na(stp_name)) %>%
group_by(stp_name, vtmnm) %>%
summarise(volume_ddd = sum(volume_ddd)) %>%
group_by(stp_name) %>%
mutate(prop_use = volume_ddd / sum(volume_ddd),
pos = cumsum(volume_ddd) - volume_ddd/2,
total = sum(volume_ddd)) %>%
ungroup() %>%
mutate(rank = dense_rank(-total),
stp_name = fct_reorder(stp_name, -rank))
temp_ggplot <- df_scmd_asthma_ddd_map_temp %>%
filter(rank <= 20) %>%
ggplot(aes(stp_name)) +
geom_bar(aes(y = volume_ddd,
fill = vtmnm,
text = paste0("<b>STP:</b> ", stp_name, "<br>",
"<b>Total volume in ddd:</b> ", round(total, 0), "<br>",
# "<b>Medication:</b> ", vtmnm, "<br>",
"<b>", vtmnm , " volume in ddd:</b> ", round(volume_ddd, 0), " (", scales::percent(prop_use, accuracy = 0.1), ")"
)
),
# text = paste0("<b>STP:</b> ", stp_name, "<br>",
# "<b>Volume:</b> ", round(total, 0), "<br>",
# "<b>Proportion:</b> ", scales::percent(prop_use, accuracy = 0.1), "<br>",
# "<b>Medication:</b> ", vtmnm)),
stat='identity',
# position = position_dodge()
) +
scale_fill_viridis_d() +
labs(subtitle = paste0("From: ", min(df_scmd_asthma_ddd_map$year_month), " to ", max(df_scmd_asthma_ddd_map$year_month)),
x = NULL,
y = "Defined Daily Dose",
fill = NULL) +
scale_y_continuous(labels = scales::comma) +
theme(text = element_text(size = 12),
legend.position="bottom") +
coord_flip()
# temp_ggplot
plotly::ggplotly(temp_ggplot,
tooltip = "text") %>%
plotly::config(displayModeBar = FALSE) %>%
layout(legend = list(orientation = "h", x = -0.5, y =-.15))Figure. Total prescribing of asthma biologics for 20 STPs with the largest volume across all selected medications.
green_pal <- function(x) rgb(colorRamp(c("#effbf7", "#31cf96"))(x), maxColorValue = 255)
df_scmd_asthma_ddd_map_temp_temp <- df_scmd_asthma_ddd_map_temp %>%
mutate(prop_use = round(prop_use, 3)) %>%
pivot_wider(id_cols = c(stp_name, total, rank),
names_from = vtmnm,
values_from = c(volume_ddd, prop_use),
values_fill = 0) %>%
select(stp_name, starts_with("prop"), total)
data_tab <- df_scmd_asthma_ddd_map_temp_temp %>%
mutate(prop_scale = seq(from = 0, to = 1, length = nrow(df_scmd_asthma_ddd_map_temp_temp)))
reactable(data_tab,
filterable = TRUE,
# searchable = TRUE,
columns = list(
stp_name = colDef(name = "STP",
minWidth = 300),
prop_use_Benralizumab = colDef(name = "Benralizumab",
format = colFormat(percent = TRUE, digits = 1),
minWidth = 100,
style = function(value) {
normalized <- (value - min(data_tab$prop_scale)) / (max(data_tab$prop_scale) - min(data_tab$prop_scale))
color <- green_pal(normalized)
list(background = color)
}),
prop_use_Mepolizumab = colDef(name = "Mepolizumab",
format = colFormat(percent = TRUE, digits = 1),
minWidth = 100,
style = function(value) {
normalized <- (value - min(data_tab$prop_scale)) / (max(data_tab$prop_scale) - min(data_tab$prop_scale))
color <- green_pal(normalized)
list(background = color)
}),
prop_use_Omalizumab = colDef(name = "Omalizumab",
format = colFormat(percent = TRUE, digits = 1),
minWidth = 100,
style = function(value) {
normalized <- (value - min(data_tab$prop_scale)) / (max(data_tab$prop_scale) - min(data_tab$prop_scale))
color <- green_pal(normalized)
list(background = color)
}),
prop_use_Reslizumab = colDef(name = "Reslizumab",
format = colFormat(percent = TRUE, digits = 1),
minWidth = 100,
style = function(value) {
normalized <- (value - min(data_tab$prop_scale)) / (max(data_tab$prop_scale) - min(data_tab$prop_scale))
color <- green_pal(normalized)
list(background = color)
}),
total = colDef(name = "Total (ddd)",
minWidth = 65,
format = colFormat(digits = 0)),
prop_scale = reactable::colDef(show = FALSE)
),
style = list(fontSize = "12px"),
highlight = TRUE)